clearvars
clc

% Produce table 1 by loading different calibrations and storing results

Results = NaN(13,4);

for ical=1:4

    name = strcat('Results',num2str(ical));
    load(name)
    
    % Construct employment in levels in the counterfactual economy
    
    LDCfa = LHCfa .* LDBas;
    LLBas = lLBas;
    LLBas(:,2:end,2:end) = repmat(lLBas(:,2:end,1),1,1,T-1).*cumprod(LDBas,3);
    LLCfa = lLCfa;
    LLCfa(:,2:end,2:end) = repmat(lLCfa(:,2:end,1),1,1,T-1).*cumprod(LDCfa,3);
    
    % Shares, welfare, and bartik
    
    lShaReg            = lLBas(:,:,2)./repmat(sum(lLBas(:,:,2),2),1,Splus);
    AggWelf            = sum(WELF .* lShaReg,2)*100;
    AggWelfU           = AggWelf(1:M);
    bartikjose         = xlsread('Inputs/InputData.xlsx','EXP','B2:B51');
    barr               = bartikjose*2.63/mean(bartikjose);
    
    % Compute baseline things
    
    PopulationB        = squeeze(sum(lLBas(1:M,:,:),2));
    TotalLaborSupplyB  = squeeze(sum(lLBas(1:M,2:end,:),2));
    TotalEmploymentB   = squeeze(sum(LLBas(1:M,2:end,:),2));
    MoPB               = squeeze(sum(LLBas(1:M,2:13,:),2))./PopulationB;
    NoPB               = squeeze(sum(LLBas(1:M,14:15,:),2))./PopulationB;
    UoPB               = (TotalLaborSupplyB-TotalEmploymentB)./PopulationB;
    loPB               = TotalLaborSupplyB./PopulationB;
    EoPB               = TotalEmploymentB./PopulationB;
    WagebillMB         = squeeze(sum(YLBas(1:M,1:12,:),2));
    LaborMB            = squeeze(sum(LLBas(1:M,2:13,:),2));
    WagesMB            = WagebillMB./LaborMB;
    WagebillNB         = squeeze(sum(YLBas(1:M,13:end,:),2));
    LaborNB            = squeeze(sum(LLBas(1:M,14:end,:),2));
    WagesNB            = WagebillNB./LaborNB;
    MoPchanbas         = (mean(MoPB(:,7:9),2)-MoPB(:,1))*10/7*100;
    NoPchanbas         = (mean(NoPB(:,7:9),2)-NoPB(:,1))*10/7*100;
    UoPchanbas         = (mean(UoPB(:,7:9),2)-UoPB(:,1))*10/7*100;
    loPchanbas         = (mean(loPB(:,7:9),2)-loPB(:,1))*10/7*100;
    EoPchanbas         = (mean(EoPB(:,7:9),2)-EoPB(:,1))*10/7*100;
    PoPchanbas         = (PopulationB(:,8)-PopulationB(:,1))./...
                          PopulationB(:,1)*10/7*100;
    MWchanbas          = (WagesMB(:,8)./WagesMB(:,1)-1)*100*10/7;
    NWchanbas          = (WagesNB(:,8)./WagesNB(:,1)-1)*100*10/7;
    
    % Compute counterfactual things
    
    PopulationC        = squeeze(sum(lLCfa(1:M,:,:),2));
    TotalLaborSupplyC  = squeeze(sum(lLCfa(1:M,2:end,:),2));
    TotalEmploymentC   = squeeze(sum(LLCfa(1:M,2:end,:),2));
    MoPC               = squeeze(sum(LLCfa(1:M,2:13,:),2))./PopulationC;
    NoPC               = squeeze(sum(LLCfa(1:M,14:15,:),2))./PopulationC;
    UoPC               = (TotalLaborSupplyC-TotalEmploymentC)./PopulationC;
    loPC               = TotalLaborSupplyC./PopulationC;
    EoPC               = TotalEmploymentC./PopulationC;
    WagebillMC         = squeeze(sum(YLCfa(1:M,1:12,:),2));
    LaborMC            = squeeze(sum(LLCfa(1:M,2:13,:),2));
    WagesMC            = WagebillMC./LaborMC;
    WagebillNC         = squeeze(sum(YLCfa(1:M,13:end,:),2));
    LaborNC            = squeeze(sum(LLCfa(1:M,14:end,:),2));
    WagesNC            = WagebillNC./LaborNC;
    MoPchancfa         = (mean(MoPC(:,7:9),2)-MoPC(:,1))*10/7*100;
    NoPchancfa         = (mean(NoPC(:,7:9),2)-NoPC(:,1))*10/7*100;
    UoPchancfa         = (mean(UoPC(:,7:9),2)-UoPC(:,1))*10/7*100;
    loPchancfa         = (mean(loPC(:,7:9),2)-loPC(:,1))*10/7*100;
    EoPchancfa         = (mean(EoPC(:,7:9),2)-EoPC(:,1))*10/7*100;
    PoPchancfa         = (PopulationC(:,8)-PopulationC(:,1))./...
                          PopulationC(:,1)*10/7*100;
    MWchancfa          = (WagesMC(:,8)./WagesMC(:,1)-1)*100*10/7;
    NWchancfa          = (WagesNC(:,8)./WagesNC(:,1)-1)*100*10/7;
    
    % Obtain differences in responses between baseline and counterfactual
    
    MoPchacha          = MoPchancfa-MoPchanbas;
    NoPchacha          = NoPchancfa-NoPchanbas;
    UoPchacha          = UoPchancfa-UoPchanbas;
    loPchacha          = loPchancfa-loPchanbas;
    PoPchacha          = PoPchancfa-PoPchanbas;
    MWchacha           = MWchancfa-MWchanbas;
    NWchacha           = NWchancfa-NWchanbas;
    X                  = [ones(M,1),barr];
    Y                  = [MoPchacha,NoPchacha,UoPchacha,-loPchacha,MWchacha,...
                          NWchacha,PoPchacha,AggWelfU];
    coeff              = (X'*X)^(-1)*(X'*Y);
    betacoeffchacha    = coeff(2,:);
    
    if lecs==7
        Targettim      = [6800.53  ; 11113.98 ; 29459.96 ; 37076.37 ; ...
                          34354.78 ; 35868.28 ; 46829.88];
        Targetsec      = [7322.697 ; 31151.50 ; 2849.315 ; -260.125 ; ...
                          11683.57 ; 7511.861 ; 1800.746 ; 17282.76 ; ...
                          18738.14 ; 100334.3 ; 7663.853 ; 6821.200];
    elseif lecs==11
        Targettim      = [5104.412 ; 8592.531 ; 23428.24 ; 29587.35 ; ...
                          27386.50 ; 28610.41 ; 37474.64 ; 30501.00 ; ...
                          -25794.0 ; 52228.55 ; 57388.52];
        Targetsec      = [11627.39 ; 42766.07 ; 4187.986 ; 2528.457 ; ...
                          21558.14 ; 13690.39 ; 3778.203 ; 24939.91 ; ...
                          29002.93 ; 136960.6 ; 15674.90 ; 6127.482];
    else
        disp('Wrong lenght of time china shock')
    end
    
    Targetsec          = Targetsec.*(sum(Targettim)./sum(Targetsec));
    
    % Obtain residuals
    
    F1                 = (betacoeffchacha(3)/0.221-1)*100;
    F2                 = (betacoeffchacha(4)/0.553-1)*100;
    F3                 = (-betacoeffchacha(7)/0.050-1)*100;
    F4                 = (changetimcha./Targettim-1)*100;
    F5                 = (changeseccha(1:12)./Targetsec-1)*100;
    
    % Population share by states in USA
    
    PopSta             = sum(lLBas(1:50,:,2),2);
    PopSha             = PopSta./sum(PopSta);
    MeanWelfShaWeighed = PopSha'*AggWelfU;

    % Store results

    Results(1,ical)    = betacoeffchacha(3);
    Results(2,ical)    = betacoeffchacha(4);
    Results(3,ical)    = betacoeffchacha(1);
    Results(4,ical)    = betacoeffchacha(2);
    Results(5,ical)    = betacoeffchacha(7);
    Results(6,ical)    = betacoeffchacha(5);
    Results(7,ical)    = betacoeffchacha(6);
    Results(8,ical)    = betacoeffchacha(8);
    Results(9,ical)    = MeanWelfShaWeighed;
    Results(11,ical)   = nu;
    if ical~=3
    Results(12,ical)   = ka;
    end
    Results(13,ical)   = delta;

    % Compute welfare gain without DNWR

    name = strcat('Results',num2str(caltype),'Flex');
    load(name)
    
    lShaReg            = lLBas(:,:,2)./repmat(sum(lLBas(:,:,2),2),1,Splus);
    AggWelf            = sum(WELF .* lShaReg,2)*100;
    AggWelfU           = AggWelf(1:M);   
    PopSta             = sum(lLBas(1:50,:,2),2);
    PopSha             = PopSta./sum(PopSta);
    Results(10,ical)   = PopSha'*AggWelfU;

end

disp(compose('%.3f', Results))
